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Abstract 

I The domain wall formulation of lattice fermions is expected to support ac- 

' curate chiral symmetry, even at finite lattice spacing. Here we attempt to 

■ use this new fermion formulation to simulate two-flavor, flnite temperature 

. QCD near the chiral phase transition. In this initial study, a variety of quark 

, masses, domain wall heights and domain wall separations are explored us- 

I ing an 8 X 4 lattice. Both the expectation value of the Wilson line and the 

chiral condensate show the temperature dependence expected for the QCD 
phase transition. Further, the desired chiral properties are seen for the chiral 
. condensate, suggesting that the domain wall fermion formulation may be an 

eff^ective approach for the numerical study of QCD at flnite temperature. 



X 

S ■ I. INTRODUCTION 

Many of the properties of low energy QCD are a direct consequence of the breaking of 
chiral symmetry by the QCD vacuum. It is expected that this spontaneous chiral symmetry 
breaking will disappear as the temperature is increased. Both the nature of this symme- 
try restoration (abrupt phase transition or continuous cross-over) and the character of the 
high-temperature quark-gluon plasma phase are active areas of both theoretical |l|J^ and 
experimental research 

An especially promising approach to the theoretical study of equilibrium properties of 
both the QCD phase transition and the high-temperature plasma phase is direct numerical 
simulation of the Feynman path integral using the methods of lattice gauge theory. The 
quantum partition function is written as a Euclidean path integral that can be studied 
ab initio using the discrete, lattice formulation of Wilson 0. While the local color gauge 
symmetry of the theory remains exact at any lattice spacing in Wilson's formulation, much 
of the theory's flavor symmetry, and especially its chiral component, is explicitly broken. 



1 



This difficulty in representing the continuum flavor symmetries in a lattice fermion for- 
mulation is a serious problem that has persisted for more than two decades. When the 
fermion action is naively discretized the low-energy fermionic degrees of freedom increase 
by a factor of 2^. This well-known "doubling" problem can only be remedied by methods 
that explicitly break the chiral flavor symmetries for flnite lattice spacing 0. The chiral 
symmetries are then recovered together with the Lorentz symmetry as the lattice spacing is 
sent to zero. The most popular of these methods are staggered 10-0 and Wilson fermions. 

Although, in principle these methods should be able to approximate the continuum 
theory in a controlled way, in practice this problem has been a formidable obstacle to lattice 
studies of the QCD phase transition. For example, the Wilson fermion formulation explicitly 
breaks all of the continuum chiral symmetries making phenomena driven by the spontaneous 
breakdown of chiral symmetry difficult to study. While staggered fermions do possess a one- 
dimensional continuous chiral symmetry at flnite lattice spacing, this formulation explicitly 
breaks the vector flavor symmetry so instead of three light Goldstone pions with mass on 
the order of the critical temperature ~ 160 MeV as found in Nature, present staggered 
simulations have masses for two of the three pions in the range 500-600 GeV, certainly too 
large. 

In addition, the subtle effects of the continuum axial anomaly which are closely connected 
with the order of the transition |10] are badly mutilated by both fermion formalisms at flnite 
lattice spacing. While the anomalous Ua{^) continuum chiral symmetry is explicitly broken 
by both formalisms, the fermion zero modes required by Atiyah-Singer index theorem are 
shifted away from zero by flnite lattice spacing effects. 

In principle, each of these difficulties can be addressed by simply working at smaller 
lattice spacing. However, present numerical methods scale very poorly as the lattice spacing 
is decreased, with the required numerical effort growing as ~ l/a^~^^ for lattice spacing a. 

Domain wall fermions (DWF) offer a new approach to the problem of including fermions 
in lattice gauge theory calculations. In this formulation, introduced by Kaplan [0,0], the 
fermionic flelds are deflned on a flve- dimensional hyper-cubic lattice using a local action. 
The flfth direction can be thought of as an extra space-time dimension or as a new internal 
flavor space. The gauge flelds are represented in the standard way in four dimensional 
space-time and are coupled to the extra fermion degrees of freedom in a diagonal fashion. 

In this paper, we use a variant of Kaplan's approach, developed by Shamir [|13|, in which 
the two four- dimensional faces orthogonal to the new flfth dimension are treated differently, 
with free boundary conditions imposed on the fermion flelds. This key ingredient allows 
a system made up of naively massive fermions to develop chiral surface states on these 
boundaries (domain walls) with the positive chirality states bound exponentially to one wall 
and the negative chirality states bound to the other. 

The two chiralities overlap only by an amount that is exponentially small in Lg, the 
number of lattice sites along the flfth direction. The resulting mixed state forms a Dirac 
4-spinor that propagates in the four-dimensional space-time with an exponentially small 
mass. Therefore, the amount of chiral symmetry breaking that is artiflcially induced by 
the regulator can be controlled by the new parameter Lg. In the — oo limit the chiral 
symmetry is exact even at flnite lattice spacing. Thus, the domain wall fermion method 
has succeeded in disentangling the chiral limit {Ls oo) and the continuum limit (a — 0). 
Furthermore, the direct computing requirement grows only linearly with Lg. 
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Here we report the first full QCD simulations using domain wall fermions in four dimen- 
sions. The properties and parameter space of domain wall fermions appropriate for a study 
of QCD thermodynamics are explored in detail. Small lattices of size 8^ x 4 were used to 
perform numerical simulations of full, two-flavor QCD at finite temperature. Preliminary 
results of this work have appeared in [|T^-[T6|. These studies have been carried out using the 
QCDSP supercomputer at Columbia |T^. Based on the work reported here, results of phys- 
ical interest have been obtained on larger lattices for a variety of observables. Preliminary 
results of these studies can be found in and will be presented in follow-on papers |]18| . 



For a detailed introduction to the subject and relevant references the reader is referred 
to Refs. |l9|-pT|, and the reviews in Refs. p2|-p5|. Earlier numerical work using domain 



wall fermions has explored the parameter space of a QCD-like, dynamical vector theory 
in two dimensions, the two flavor Schwinger model p!9|j20|] . For applications to quenched 
QCD see Refs. p6|-p^,p!5|] for applications to four- Fermi models see Ref. ||38| and for possible 
alternatives to domain wall fermion simulations see Refs. PU|-^. 

In Section H the action of the theory and a brief description of the numerical methods are 
presented. In Section some important analytical facts are outlined in order to help guide 
the numerical investigation. In Section |^ we study the chiral properties of the theory both 
below and above the chiral phase transition. Our numerical results suggest that domain 
wall fermions are able to sustain the desired chiral properties of QCD, even at finite lattice 
spacing. Both a low temperature phase where the SU (2) x SU (2) chiral symmetry is broken 
spontaneously to an SU (2) vector symmetry and a high temperature phase where the full 
SU (2) X SU{2) chiral symmetry is intact can be recognized. 

In Section ^ the dependence on the two new regulator parameters, the number of sites 
in the fifth direction Lg, and the domain wall "height" mo, is studied numerically. Finally, 
in Section |V| conclusions and outlook are presented. Appendix ^ gives the exphcit form of 
the gamma matrices used in this work while Appendix B describes the molecular dynamics 
equations of motion. Tables summarizing the numerical results are given at the end of the 
paper. 



II. HYBRID MONTE CARLO WITH DOMAIN WALL FERMIONS 

In this section the action of QCD with domain wall fermions, its implementation for 
the Hybrid Monte Carlo (HMC) algorithm, and the parameters used in the simulations are 
described. In the following, we discuss the case of two degenerate flavors implemented using 
the HMC $ algorithm |^6[. (An odd number of flavors can be simulated using the HMC R 
algorithm [0). 

Domain wall fermions can be used in numerical simulations in a fashion similar to tra- 
ditional Wilson fermions. In fact, if the fifth direction is thought of as an internal flavor 
direction then an HMC simulation with DWF is identical to a simulation of many flavors 
of Wilson fermions with a sophisticated mass matrix. We use the partition function of 
QCD with domain wall fermions proposed in PTf but with a slightly different heavy flavor 



subtraction as in ||T9| , ^ . In particular: 

Z = J [dU] J [rf^rf^] J [d^^pyd^pv]e-^ (1) 
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U^{x) is the gauge field, s) is the fermion field and ^pv{x, s) is a bosonic, Pauli-Villars 
field. The variable x specifies the coordinates in the four- dimensional space-time box with 
extent L along each of the spatial directions and extent Nt along the time direction while 
s = 0,1, ... ,Ls — I is the coordinate of the fifth direction, with Lg assumed to be even. The 
action S is given by: 

S = Sg{U) + Sp{^, ^, U) + Spvi^U^ "^pv, U) (2) 

where: 

5G = /5E(l-^ReTr[f/,]) (3) 
p "J 

is the standard plaquette action, j3 = Q/g^ and is the lattice gauge coupling. The fermion 
action for two flavors is: 

Sf = - J2 ^f{x,s)DF{x,s;x',s')^f{x',s') (4) 

x,x' ,s,s',f 

with flavor index / = 1, 2 and Dirac operator: 

Df{x, s; x', s') = 6s,s'D^^ [x, x') + D^{s, s')S,,,> (5) 



D\x,x') = - E [{1 -lf,)Uf,{x)6^+f,^^> 

+ (mo - 4)(5a;,a;' 



il+^,)Ul{x% 



X — fl,X' 



(6) 
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D {S, S') = - (1 - l5)Ss+l,s' + (1 + 75)5s-l 



2 
2 



2(5. 



(1 - 75)5l.-1,s'5s',0 + (1 + 75)'5,,o'5l.-i,, 



(7) 



Here, s and s' lie in the range < s, s' < — 1. In the above equations mo is a flve- 
dimensional mass representing the height of the domain wall in Kaplan's original language. 
In order for the doubler species to be removed in the free theory one must choose < mo < 2 
The parameter rrif explicitly mixes the two chiralities and, as a result, controls the 
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bare fermion mass of the four- dimensional effective theory. 

While the DWF Dirac operator deflned above is not hermitian, it does obey the identity 
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(8) 



with i?5 the reflection operator along the flfth direction. As a result the single-flavor Dirac 
determinant is real: detDp* = deiD^ = det 'j^R^Dp^sR^ = detDp and the two-flavor 

2 



determinant which follows from integrating out the fermions in Eq. |1], det Dp , is positive. 
The gamma matrices used in this work are given in Appendix 0. Also notice that Dp is the 



same as the Dp of P7|. 
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The Pauli-Villars action is designed to cancel the contribution of the heavy fermions in 
the large Lg limit. Normally, such heavy fermions decouple from low energy physics and can 
be safely ignored. However, in the present situation the number of heavy fermions grows 
proportional to Lg and can potentially overwhelm the effects of the fixed number of low 
energy degrees of freedom of interest. Specifically this difficulty will arise for the order of 



limits for which DWF are intended: first ^ oo followed by a — > 0. ||48H51|. 

There is some fiexibility in the definition of the Pauli-Villars action since different actions 
can easily have the same — > cxd limit. However, the choice of the Pauli-Villars action may 
affect the approach to the Lg —>■ oo limit. A slightly different action than that proposed by 
Furman and Shamir |^ is used here. This action |T9| , pO| is easier to implement numerically 
and, even for finite Ls, it exactly cancels the fermion action when ruf = 1 resulting in a pure 
gauge theory. For two fermion fiavors, the Pauli-Villars action we use is: 

Spv = J2 ^Pv{x^s)MF{x,s;x',s')\„,^=i^Pvix',s') (9) 

x,x' ,s,s' 

where Mp = D^pDp. 

The traditional HMC $ algorithm was constructed directly from the action of Eq. ^. In 
order to improve performance a standard even-odd preconditioning |^ of the Dirac operator 



Dp was employed. The even-odd preconditioning was done on the five dimensional space. All 
necessary matrix inversions were done using a standard conjugate gradient (CG) algorithm. 
As expected the even-odd preconditioning resulted in a reduction of the required number of 
conjugate gradient iterations and a consequent speed-up of a factor of approximately two. 

The only new ingredient in our HMC algorithm is the appearance of the bosonic Pauli- 
Villars fields. The probability distribution of these fields is generated with a heat bath step at 
the beginning of each HMC "trajectory": a field of Gaussian random numbers is generated 
with distribution e~^pv^^^ and from it the Pauli-Villars fields $py(x,s) are obtained by 
^PV = [Dp{mf = l)]~^ripv using the CG algorithm. 

Since the Pauli-Villars action in Eq. ^ is polynomial in the domain wall operator Dp, 
its gradient with respect to the gauge fields, needed to evolve the gauge degrees of freedom, 
can be computed without performing any Dirac inversions. This contrasts favorably with 
the fermion contribution to the gauge force which requires one inversion per molecular 
dynamics step. As a result, the relative computational cost involved in calculating the 
Pauli-Villars force is negligible. Furthermore, because the Pauli-Villars fields are bosonic 
their molecular dynamics force term enters with an opposite sign that of the fermion force, 
resulting in a large, approximate cancellation. Because of this cancelation the HMC force 
term is approximately independent of Lg and it is not necessary to decrease the HMC step 
size as Lg is increased. 

In the approach described above the presence of the Pauli-Villars fields increases the 
memory requirement. However, it should be noted that there is an alternative approach 
that does not involve Pauli-Villars fields. To see this consider the result after integra- 
tion over both the Pauli-Villars and fermion fields. It is det Mp{mf) / det Mp{mf = 1) = 
det [Mp{mf)/Mp{mf = 1)]. Therefore, one could simulate the same action without Pauli- 
Villars fields by simply using as the fermion matrix Mp{mf)/Mp{mf = 1). Inversion of this 
matrix will involve inversion of Mp{mf) using the CG algorithm as in the previous method 
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while the final result would have to be multiplied by the matrix Mplmf = 1). If, for exam- 
ple, the CG algorithm required 100 iterations to converge, this extra matrix multiplication 
will increase the computing cost by only 1%. The only disadvantage of this approach is that 
the equations of motion become slightly more complicated. 

Since this work is the first to implement DWF in dynamical QCD the approach with 
Pauli-Villars fields was adopted for simplicity and because it has been proven reliable in 
numerical simulations of the Schwinger model [0,^]. For the convenience of the reader the 
molecular dynamics equations of motion with Pauli-Villars fields and an even-odd precon- 
ditioned DWF Dirac operator are given in Appendix O. 



Fermionic Green's functions were computed using the method described in Ref. p7 
Standard fermion fields in the four-dimensional space-time are constructed from the five- 
dimensional fermion fields using the projection prescription: 

ijix) = Pl*(x, 0) + Pr^(x, Ls - 1) 

= ^(x, - 1)Pl + ^(x, 0)Pr (10) 

where Pr/l = |(1 i 7^)- This somewhat arbitrary choice defines operators which should 
have a large overlap with the physical low energy fermion modes bound to the s = and 
s = Lg — 1 walls. The right- and left-handed components found on opposite walls are 
combined to assemble the desired physical 4-spinors. 

Since these are the first simulations of DWF in dynamical QCD there are no previous 
results that would allow an independent check of the methods and code. Tests using the chi- 



ral condensate from the free field analytical results of |I9|J20|] were done in order to check the 
Dirac operator and inverter. The subtraction of Pauli-Villars fields was tested by performing 
simulations with rrif = 1 and comparing with equivalent results from quenched simulations. 
Finally, two flavor dynamical simulations were done on 2^ lattices and the results were com- 
pared with simulations using the overlap formalism P5|-^ relevant for the DWF action 
^ for the same parameters. In particular for (3 = 5.6, ruf = 0.1, mo = 0.9 the overlap 
simulation gave (ipip) = 1.672(2) x 10"'^ and average plaquette (plaq) = 5.765(79) x 10~^ 
while the DWF simulation with Lg = 18 gave (ipip) = 1.653(33) x 10"'^ and average plaquette 
(plaq) = 5.841(47) x 10"^. 

All numerical results in this work were obtained from lattices of size L = 8, iV^ = 4 with 
periodic spatial boundary conditions and anti-periodic temporal boundary conditions. The 
fifth direction was set to various values in the range [8,40], the domain wall height mg was 
varied in the range [1.15,2.4], the fermion mass was varied in the range [0.02,0.18] and f3 
was varied in the range [4.65,5.95]. The molecular dynamics trajectory length was set to 
r = 0.5 and the step size 6t was set to various values in the range [0.0078, 0.02] depending 
on the values of the other parameters. The CG stopping condition which is defined as the 
ratio of the norm of the residual vector over the norm of the source was set to 10~^. This 
resulted in an average number of CG iterations ranging between 50 and 400 depending on 
the values of the other parameters. 

The initial configuration was generally chosen to be in the phase opposite to that expected 
for the input parameters creating a very visible thermalization process in which the system 
should be seen to evolve into the correct phase. Typically 100 — 400 trajectories were needed 
to thermalize the lattice. The chiral condensate and Wilson line were measured in every 
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sweep. The chiral condensate was measured using a standard "one-hit" stochastic estimator 
of the trace of Dp^ with spin and s coordinates restricted according to Eq. |l^. Specifically 
we evaluated the quantities: 

(l^l) = TT^|Etr[ n U,]\ (11) 

X l£L{x) 

^^"^^ = 12Z^{*'[^' = 0|1/Z}^|s = U- 1)Pr] 

+tr[(s = L, -1|1/D^|s = 0)Pl]}. (12) 

Here Ui identifies the SU (3) gauge matrix corresponding to the link I and the ordered product 
is taken for all links in the time-like line L{x) with spatial coordinate x. The somewhat 
unconventional normalization in Eq. |12| was used in our previous work and determines a 
spin and color average which for very large mass mj approaches l/ruf. (Note, here Dp is 
the single-flavor Dirac operator defined in Eq. ^.) 



III. ANALYTICAL CONSIDERATIONS 

In this section we summarize some of the analytically determined properties of domain 
wall fermions. These help guide our numerical investigations, which are done for finite and 
non-zero values for the three parameters of domain wall fermions, Lg, rriQ, and rrif, as well 
as at finite bare coupling go. 



A. Ls dependence 

For numerical simulations, the existence of the chiral limit for domain wall fermions and 
the rate of approach to it are of primary importance. The computational requirements for 
domain wall fermions grow as one power of Lg from the simple increase in the number of 
operations. An additional slight increase in computational cost for larger Lg comes from the 
decrease in the total quark mass due to smaller mixing between the chiral surface states, 
until the quark mass is dominated by the input rrif. 

The axial Ward-Takahashi identities for domain wall fermions are the same as the contin- 
uum, except for an additional term which comes from the mixing of the left- and right-handed 
light surface states at the midpoint of the fifth dimension, Ls/2. At any lattice spacing this 
additional term vanishes as — >■ oo for non-singlet axial symmetries P7|-pl|] . For the singlet 
axial symmetry, this extra term generates the axial anomaly. At strong coupling, the axial 
currents are conserved for Lg ^ oo but, since the doubler fermions may enter the spectrum. 



these currents may not have the physical significance of axial currents . 

For free domain wall fermions, the rate of approach to the chiral limit can be calculated. 
At finite L<, the mixing of the chiral components is reflected in the fermion mass rries- For 



the one flavor theory this effective mass is ||20|| 



rrics = ''^0(2 — mo) mj + (1 — rno)^" < mo < 2, (13) 
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meff has two pieces: one is proportional to the bare mass mj and the other expresses the 
residual mixing between the chiral modes bound to the domain walls. Since each bound chiral 
state decays exponentially with the distance from its wall, the residual mixing between them 
vanishes exponentially with Lg, with a decay constant of — In |1 — mo|. Notice that when 
Lg oo, mo becomes an irrelevant parameter, provided it stays in the range (0,2). 

In the free theory, one also finds that fermion states with non-zero four-momentum decay 
more slowly with the distance from the wall than do zero momentum states. The decay 
is controlled by the four-momentum and the value for rriQ. Since the lattice momentum 
Pl — P^^i where a is the lattice spacing, the slower decay for modes with non-zero four- 
momentum is an 0{a?) effect which should vanish in the continuum limit. In addition, for a 
given mo, there is a critical four-momentum above which the fermions are no longer bound 
to the wall, but instead behave like massive, five- dimensional fermions. Of course, because 
these fermions are massive, they necessarily preserve the theory's four- dimensional chiral 
symmetry since their propagation between the s = and s = Lg — \ walls is exponentially 
suppressed. 

For interacting theories, a simple expectation is for Eq. ^ to be replaced by 

meff = Z„ [m/ + c e-"^=] . (14) 
The exponential dependence is seen perturbatively [l^JS^JS^ and proven to exist non- 



perturbatively, provided the gauge fields satisfy a smoothness condition ||55| , |56| . These ana- 
lytic results support the expectation of exponential suppression of chiral symmetry breaking 
effects in the non-perturbative regime. However, this behavior may be best established by 
the sort of explicit numerical study reported here. Generally a should depend on mo, al- 
lowing one to choose an optimal value for simulations at finite Lg. While in the free theory 
mo = 1 gives e~" = 0, for the interacting theory the variable character of fermion propaga- 
tion in fluctuating background gauge fields makes decoupling the walls with a single value 
for mo unhkely, except at very weak coupling. 

Close to the continuum limit, it can be argued that this form for the effective mass, 
an input quark mass plus a residual mass mres, should enter all long-distance observables. 
However, away from the continuum limit or for quantities that cannot be obtained from a 
low energy effective QCD Lagrangian this is not necessarily the case. Therefore, different ob- 
servables may approach their Lg —>■ oo limit in different ways, depending on the momentum 
scales which enter the observable, and the corrections to the input quark mass, particularly 
at stronger couplings, may be more complicated. In a numerical investigation this has to be 
kept in mind. In this paper only the chiral condensate and pion susceptibility are considered. 
Work on larger lattices involving measurements of many fermionic operators is currently in 
progress [p!8| . 

Numerical simulations may well be the only way to determine the dependence of chiral 
symmetry breaking effects on Lg for intermediate lattice spacings (~ 1 to 3 GeV). While for 
full QCD, perturbative and non-perturbative arguments support exponential falloff with Lg, 
for quenched theories, where the lack of damping from a fermionic determinant can lead to 
configurations with unsuppressed small eigenvalues for the fermions, the large Lg behavior is 
even more in need of determination through simulations. Some results from quenched QCD 



simulations have been discussed in Refs. E6hB7,15 
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B. rriQ dependence 



For free domain wall fermions the number of light flavors is controlled by the value of mo 
r^ . In particular ttlq < corresponds to zero light flavors, < itlq < 2 to one, 2 < mo < 4 



to four, and 4 < mo < 6 to six light flavors. The theory is symmetric under mo 10 — mo- 
For the interacting theory the values of mo which distinguish between different numbers of 
flavors are changed. Light fermions first appear for mo > 0, the one to four flavor transition 
occurs for mo > 2, etc. and the theory is still symmetric about mo — 10 — mo. This 
is expected perturbatively and seen numerically [^|T5tp!6| . There is also some numerical 



evidence that the transition between different numbers of flavors is smooth and spread out 
over a small region of mo [^. For the interacting theory, keeping mo < 2 guarantees that 
a theory with not more than one flavor is being studied. 

While mo is an irrelevant parameter for ^ cxd, it is very important for simulations, not 
only in controlling the approach to the chiral limit and the flavor content of the theory, but 
also for insuring that light fermions with an average momentum given by the temperature 
are still bound to the walls. For the free theory, the range of four-momenta carried by states 
that are bound to the walls increases as mo increases from zero, as do the corresponding 
Dirac eigenvalues. As mo approaches one, the largest Dirac eigenvalues of these "bound" 
states become farther off-shell, with values ~ 1/a. As mo increases above 1, the number of 
these off-shell states continues to grow but rather than their eigenvalues increasing, instead 
their degeneracy increases beyond what would be seen for the large momentum states of 
a free theory. As mo increases further and approaches 2, some of these excess, degenerate 
states become more nearly on-shell until for mo > 2 one has the low-lying Dirac eigenvalues 
of a free, four-flavor theory. Thus, in the free case a choice of mo midway between and 2 is 
best, giving the largest phase space for physical states bound to the walls, without adding 
additional flavors. Using this behavior as a guide for the interacting case, one expects that 
choosing mo midway between the value where a single light fermion is bound to the walls 
and four light fermions are bound allows the largest range of four-momentum for a single 
flavor of light quark bound to the walls. 



C. Topology 

An important property of the domain wall fermion Dirac operator is the presence of 



exact zero modes in the Lg = oo limit, as can be seen from the overlap formalism ||48H51 



These zero modes are related to the topological charge of the gauge field and as a result an 



approximate form of the index theorem is present on the lattice [57|. Studies on semiclassical 



configurations show the presence of modes which are very close to zero modes even at finite 



[21] and as a result make lattice studies of anomalous symmetry breaking possible [14-16 
During simulations, field configurations of different winding number should show zero 
mode effects in fermionic observables. The efficiency with which the hybrid Monte Carlo 
can move the system between sectors of different winding is an important question, as are 
the long correlations along the fifth direction which develop for gauge field configurations 
where the topology is changing. These issues have been studied in numerical simulations 



of the dynamical Schwinger model [BO] where the hybrid Monte Carlo performed well and 
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topology changing occurred. For this exploratory study of full QCD thermodynamics, the 
input quark masses are not small, so the effects of topology should not be particularly large. 



IV. THE FINITE TEMPERATURE QCD PHASE TRANSITION 

The previous sections have described the domain wall fermion formulation and important 
questions about it that need to be investigated numerically. Here we report on simulations 
of full QCD at finite temperature with domain wall fermions on 8^ x 4 lattices. Studying this 
system allows us to investigate domain wall fermions for full QCD and look for the presence 
of chirally broken and symmetric phases. The small volume makes scanning over many 
values for Lg, mo, rrif and go possible, laying the foundation for more realistic simulations 
on larger volume. 

Since the finite temperature transition of QCD is controlled by the chiral symmetries 
of the theory (for light quarks), using domain wall fermions to preserve the full global 
symmetries of the continuum should remove one systematic lattice error that is difficult to 
control. However, finite temperature simulations are generally only possible on relatively 
coarse lattices (a~^ ~ 700 MeV for a lattice with A^^ = 4), where analytic results about 
domain wall fermions are lacking. The light chiral modes of domain wall fermions at weak 
coupling must exist at ~ 700 MeV, in the full non-perturbative gauge field backgrounds, 
for thermodynamic simulations to be possible. If it is found that chiral modes exist on coarse 
lattices, the size of the mj-es{Ls,P) and its dependence on Lg and (3 must be investigated. 
(As already mentioned, m,.es is only a sensible quantity for low-energy observables and it 
must be demonstrated that various determinations of it are consistent. In this section we 
refer to m^es, without specifying precisely how it may be determined, as a generic indicator 
of the mixing between the chiral modes.) 



A. Locating the transition 

Locating the phase transition in full QCD requires scanning values for four parameters 
(mo, Ls, nif, and /5). Without any knowledge of the location of the transition, or if it even 
exists for domain wall fermions, choosing parameters for initial simulations is difficult. For 
staggered fermions, the critical coupling for the finite temperature phase transition for 2 
fiavors on an iVf = 4 lattice is pc = 5.265 for m = 0.01 and Pc = 5.291 for m = 0.025 



Since staggered and domain wall fermions both have their chiral limit at zero quark mass, 
the light quarks have the largest effect in the location of (3c and both theories have the same 
number of light fiavors, we used the staggered values as a rough guide. 

Our first simulations with domain wall fermions were done at /? = 5.0 and (3 = 5.4, with 
the hope that these would be above and below the transition region. Lg = 8 and rrif = 0.1 
were chosen to keep the computational difficulty modest. We worked with mo = 1.65, since 
for quenched simulations this choice gave a reasonable falloff between the walls at (3 = 5.7 
and for quenched QCD, /3 = 5.7 is close to Pc for an Nt = 4 lattice. 

Although with this choice of mo, the /? range being examined (5.0 < P < 5.4) lies below 
the chiral transition, we describe this point first since it demonstrates our very first efforts 
in charting this parameter space and the difficulties we encountered. The evolution of (ipip) 
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for ruf = 0.1 and f3 = 5.0 is shown in the upper panel of Figure |l| and the lower panel 
is for f3 = 5.4. The hybrid Monte Carlo was run with a step size of 6t = 0.025 and 20 
steps per trajectory, giving an acceptance of 66% for (3 = 5.0 and 70% for P = 5.4. The 
evolution appears quite generic and the simulation presented no difficulty to the hybrid 
Monte Carlo. For j3 = 5.0 the Wilson line expectation value was 0.0223(15) and for (3 = 5.4 
it was 0.0466(41). Both these values are small and indicate that both j3 values correspond 
to the confined phase. 

The chiral condensate (ipip) was also measured for a variety of valence masses. In 
quenched QCD at zero temperature, extrapolations of (il'ip) to rrif = using quark masses 
from rrif = 0.02 to = 0.1 were used to see that chiral modes existed for a particular mo 
p8| . The limit {il)ip){mf 0) could only be non-zero if light chiral modes were present, 
provided Ls is large enough that the residual mixing is unimportant. (For the current finite 
temperature case, {il)ip){mf — >■ 0) can be zero either from the absence of chiral modes or 
because the system is in the symmetry-restored phase.) 

Figure |^ shows that {ipip) extrapolates to a non-zero value for both (3 = 5.0 and (3 = 5.4 
and this value is not very sensitive to (3. The values for the Wilson line indicate both f3 
values are in the confined phase, so the {ipil') results show that light chiral modes are present 
with an unknown residual mixing. The insensitivity to j3 is an interesting feature. 

Next, instead of scanning larger values of /3, we decided to change ttlq from 1.65 to 1.9, 
keeping all other parameters identical. (This reflects our initial search path in parameter 
space and does not imply the absence of a transition at mo = 1.65 and Lg = 8.) The 
acceptance is 59% for /3 = 5.0 and 71% for (3 = 5.4 The Wilson line for (3 = 5.0 is 0.030(2), 
while for (3 = 5.4 it is 0.202(5), indicating that f3 = 5.4 is likely deconfined. The evolutions 
show a very different behavior for the condensate evaluated at the dynamical quark mass. 
The value at f3 = 5.0 has increased, part of which likely reflects the change with mo in the 
overlap between the flve-dimensional light modes and the surfaces at s = and Lg — 1. The 
f3 = 5.4 values are much smaller, consistent with the deconflned phase. 

Figure ^ shows the valence quark extrapolation. The small value of (■?/''?/') (m/ —>■ 0) 
suggests the restoration of chiral symmetry. Of course, there is a possibility that this small 
value might instead be caused by the loss of chiral modes. However, this is unlikely because 
we have seen that chiral modes do exist for (3 = 5.0 and one expects that at the weaker 
f3 = 5.4 coupling these chiral modes should be even more numerous. Therefore, we have 
preliminary evidence for two phases of full QCD with dynamical domain wall fermions. 

To solidify the evidence for two different phases of QCD with domain wall fermions fur- 
ther simulations for Lg = 8 and mo = 1.9 were done with dynamical quark masses of 0.14 and 
0.18. These points are shown in Figure ^. The dashed line is the flt to the quenched extrap- 
olation shown in Figure |. There is not a large difference between the two extrapolations, 
although both full QCD extrapolations fall below the quenched extrapolations, indicating 
some suppression of small eigenvalues through the presence of the fermion determinant. In 
the next section, we study the dynamical mass extrapolation of (ipip) for larger values of Lg 
to see if the non-zero value for {il)il)){mf 0) decreases with increasing Lg. 

Additional simulations with mj = 0.1, Lg = 8 and mo = 1.9 were done for f3 = 5.2, 
5.3 and 5.45, which produced the data for {ipip) and the Wilson line shown in Figure |^. 
Crossover behavior is seen for both observables further supporting the identiflcation of both 
a chirally broken and a chirally restored phase. These simulations are at a small value of 
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Ls, so the contribution of mres to the effective quark mass may be large. Since mi.es(/3, Ls) is 
hkely varying across the transition region, due to the change in /?, the shape of the curves 
is expected to reflect this varying effective quark mass. 



B. Ls dependence in the two phases 

With this evidence for two phases, we turned to exploring the Lg dependence in each 
phase. For the confined phase, we chose (3 = 5.2 to be at weaker coupling while still in this 
phase and in the deconfined phase we chose f3 = 5.45, to be farther from the transition. 
Keeping rriQ = 1.90, simulations were done for many values of Lg and the dynamical quark 
mass, TUf. Table |I| gives the parameters for f3 = 5.2 and Table || gives them for /3 = 5.45. 
A plot of the evolution of (ipip) for f3 = 5.20 and 5.45 is shown in Figure ^ for mj = 0.02 
and Ls = 16. With a step size of 5r = 1/64 the acceptance was 90%. Once again there is 
no evidence for difficulty in the hybrid Monte Carlo evolution of this system. 

Figure ^ shows results for (ipip) at /5 = 5.2 plotted versus rrif for L<, = 8 and 16. The 
dashed lines are linear fits to the lowest three values for rrif while the solid lines are quadratic 
fits to all values oirrif. The fits for Ls = 8 are 

{^pi;) = 0.0117(2) + 0.095(2)m/ (15) 
{ipi;) = 0.0112(3) + 0.114(5)m/ - 0.15(2)mJ (16) 

with A^dof = 1 and 2 and x^/^doi = 3.7 and 0.4, respectively. The fits for Ls = 16 are 

{ipi;) = 0.0082(1) + 0.089(2)m/ (17) 
{^P^) = 0.0080(2) + 0.099(3)m/ - 0.08(l)mJ (18) 

with A^'dof = 1 and 2 and x^/Adof = 0.03 and 0.5, respectively. The results shows a strong 
Ls dependence to which we now turn. 

Figure |^ shows (ipip) for /? = 5.2 plotted versus Ls for a variety of values of rrif. The 
curves are fits to the form cq + Ci exp(— aL^) for = 8 to 40. The fit parameters are 

{^p1p) = 0.01527(4) + 0.0188(8) exp(-0.149(5)L,) mj = 0.1 (19) 
(xp^) = 0.00779(8) + 0.014(1) exp(-0.116(8)L,) mj = 0.02 (20) 
(V^V') = 0.0059(1) + 0.014(1) exp(-0.11(l)L,) mj 0.0 (21) 

All fits have Adof = 4 and give x^/^dof = 5.1, 5.6 and 6.6, respectively. The rrij 
points are first found by extrapolating to m/ = at fixed Ls and then fitting these values 
versus Ls- Although the values for are somewhat large, the data is well fit by a function 
with exponential dependence on Ls. (Note these somewhat large values can be caused 
by underestimates of the errors which may result if our Monte Carlo evolutions are not 
sufficiently long to allow proper control the long-time autocorrelations.) 



Similar results have been obtained for (3 = 5.45. Figure shows the results for {i/jiIj) 
for f3 = 5.45 for Ls = 8 and 16. {Ls = 24 and 32 results are tabulated below.) Again, 
the dashed lines are linear fits to the lowest three values for rrif while the solid lines are 
quadratic fits to all values of rrif. The fits for Lg = 8 are 
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= 0.00227(7) + 0.095(l)m/ (22) 

{^p'^p) = 0.00219(9) + 0.099(2)m/ - 0.037(9)mJ (23) 

with A^dof = 1 and 2 and x^/^doi = 0.6 and 0.1, respectively. The fits for = 16 are 

{tptp) = 0.00039(8) + 0.100(2)m/ (24) 

i^ptp) = 0.00040(6) + 0.100(3)m/ - 0.01(2)m5 (25) 

with Nfiof = 1 and 2 and x^/^doi = 0.09 and 0.02, respectively. Linear fits for the larger 
values of Ls give 

{tptp) = 0.00016(8) + 0.100(2)m/ = 24 (26) 

{ipip) = 0.00006(6) + 0.099(l)m/ = 32 (27) 

with A^dof = 1 for both Lg and x^/^dof = 0.01 and 7.1, respectively. We see that with 
increasing Lg, the extrapolated value for the condensate at mj = decreases steadily. 

Figure |rT| shows (ipip) for f3 = 5.45 plotted versus Lg for a variety of values of mj. The 
curves are fits to the form cq + Ci exp(— aL^) for = 8 to 32. The fit parameters are 



(ipij) = 0.0102(1) + 0.08(3) exp(-0.48(6)L,) mj = 0.1 (28) 

(tptp) = 0.00599(4) + 0.015(3) exp(-0.26(2)L,) ruf = 0.06 (29) 

{ipip) = 0.00213(4) + 0.025(4) exp(-0.31(2)L,) m/ = 0.02 (30) 

{ipip) = 0.00010(5) + 0.019(3) exp(-0.27(2)L,) nij 0.0 (31) 



All fits have A^dof = 3 and give x^/^dof = 0.4, 4.8, 1.1 and 0.8, respectively. Here again the 
data strongly support exponential suppression of mixing between the walls for (tptp) . 

For both the confined and deconfined cases, we see {ijjip) exponentially approaching a 
limiting value for large Lg (which is zero in the deconfined case). At the stronger coupling 
of the confined phase, the decay constant is ~ 1/10, while in the deconfined phase it is 
~ 1/4. One expects faster decay at weak coupling, but at present we do not know whether 
the different phases also play a role in the decay constant. 



C. Studying the mo dependence of the transition 

The parameter mo is relevant at finite lattice spacing, since it controls not only when 
there is a single light fermion bound to the domain walls but also the maximum momentum 
this fermion can have while still being bound. It is expected that this parameter will not 
have to be fine-tuned for domain wall fermions to work correctly, but care in choosing a 
value is necessary to get the correct number of light species and the maximum allowable 
phase space for light fermions in the thermal ensemble. 

We have studied the characteristics of the transition region by choosing rrif = 0.1, Lg = 
12 and simulating for values of /3 near the phase transition for mo = 1.15, 1.4, 1.65, 1.8, 1.9, 
2.0, 2.15 and 2.4. Tables |III|, 0, 0, ^ |VIT|, |X| and ^ contain simulation parameters 



and results. For parameters where a deconfined thermal state was expected, the initial 
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lattice was disordered, while an initial ordered lattice was used where a confined state was 
expected. 

Figure |1^ shows the expectation value of the magnitude of the Wilson line {\W\) for these 
runs. A rapid crossover is seen for all values of rriQ. The lines are the result of fitting the 
four points nearest the transition (five points where we have a point close to the transition) 
to the function 

/(x) = co(ci + tanh[c2(x-/3e)]). (32) 

This is a phenomenologically useful form for determining the point of maximum slope for 
the Wilson line. The points far from the transition are not included in these fits, since this 
phenomenological function poorly represents the data there. 

Figure |l^ shows similar results for {ijjijj) with the lines being a fit to Eq. For mo =1.15 
and 1.4, the (ipip) data do not allow even a rough determination of jSc- For small enough 
mo, the light chiral modes should not exist and we have evidence for that at mo = 1.15. 
The value for (ipip) is very small and shows little change even when the Wilson line shows 
evidence for the transition. In addition, the Wilson lines indicate the transition is very close 
to the value of 5.6925 for quenched QCD on a 24^ x 4 lattice supporting the conclusion 



that light fermion modes are not present in the simulations. The effects of the heavy modes 
are apparently quite well canceled by the Pauh-Villars fields. 

Figure gives estimates for f3c determined from the Wilson line and {ipil')- These are 



in quite reasonable agreement, particularly given the phenomenological character of their 
determination. For mo ~ 1.2, jSc is close to the quenched value and moves smoothly to 
smaller values as mg is increased. For these larger values for mo, the light quark states 
appear and the maximum momentum for a state bound to the walls should increase. These 
light states make (ipip) show crossover behavior and are required for our simulations to be 
proper studies of two-fiavor QCD. At our largest value of mo (2.4), we may be approaching 
the transition from a two fiavor theory to an eight fiavor one (recall that the domain wall 
determinant is squared in our simulations, doubling the number of fermion fiavors.) 



V. DETERMINING THE RESIDUAL MASS 



As mentioned in Section |T|, it can be expected that for long-distance physical quantities, 
the effects of mixing between the chiral wall states will result in a residual mass contribution 
to the total quark mass. This is just the statement that the dominant effect of the mixing, 
from the perspective of a low-energy effective Lagrangian, is to introduce another source 
for chiral symmetry breaking (beyond the input mj), which takes the form of the operator 
mT-esi/jiJj at low energies. For a quantity like m^, whose dependence on chiral symmetry 
breaking can be expressed as a physical parameter times the total quark mass, the quark 
mass which enters should be mj + mres- 

However, for quantities whose sensitivity to chiral symmetry breaking effects extends up 
to the cutoff scale, such an argument does not go through. The chiral condensate, {iJip) is 
such a quantity. For domain wall fermions with Lg ^ oo (or staggered fermions), expanding 
in the input quark mass in the chirally broken phase gives 
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(V'V') = Co + ciruf + 0{m'j). 



(33) 



The coefficient Ci is ultraviolet divergent in the continuum and therefore, on the lattice, gets 
large contributions from modes at the cutoff scale. For such an operator, the Lg dependence 
is not reliably represented by just making the replacement rrif rrif + mres- 

From this discussion, it is clear that although Figure ^ shows that the large Ls limit 
for (tptp) a.t rrif = 0.02 has likely been reached by Lg ~ 40, one cannot conclude that the 
value for mres has vanished. To measure mres, it is natural to look for effects in the pion 
mass, which is in turn governed by the axial Ward-Takahashi identity. This has been done 
in quenched simulations Refs. [P'^- PD| , P^| - P^ , P^ , P^ , P!^ ] , at zero temperature, but here we are 
interested in determining m^es in the confined phase at finite temperature for small volumes 
for Nf = 2 QCD. 

Our small volumes preclude taking large separations in two-point functions to completely 
isolate the pion from other states. Thus a direct measurement of the pion mass or the overlap 
of the pion with any particular source is not possible here. Instead, we use the integrated 
form for the flavor non-singlet axial Ward-Takahashi identity and try to see the contributions 
of the pion. In the zero quark mass limit on infinite volumes, the pion contributions become 
poles. Thus we can look for the effects of these precursors of the pion poles, even when they 
do not completely dominate the Ward-Takahashi identity. 

Starting from the flavor non-singlet axial Ward-Takahashi identity in |57] and summing 
over all lattice points gives 

i^Pi;) = rrifx. + AJ5. (34) 

Here ip is the four-dimensional fermion field defined by Eq. |l^ and the pseudoscalar suscep- 
tibility is (no sum on a) 

X. = ^ E ^^(^)75y^(x) ^(0)75y^(0)^ , (35) 

(The factor of l/4A^c is needed to match our normalization for (ipip).) The additional 
contribution from chiral mixing due to finite Lg is 

= ^ E (^f, (x, Lg/2) ^(0)75y ^(0)^ , (36) 

where 

r,{x, Lg/2) = ^^(x, L,/2)(l - 75)A''vl/(x, L,/2 + 1) (37) 

- ^^(x, Lg/2 + 1)(1 + 75)A"^(x, Lg/2) 

is a pseudoscalar density at the midpoint of the fifth dimension which couples left- and 
right-handed degrees of freedom. 

We have done extensive simulations for many values of Lg with f3 = 5.2, mo = 1.9 and 
m/ = 0.02 to study the consequences of the Ward-Takahashi identity. At the time of these 
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simulations, we were not measuring AJ5 explicitly. However, the other two terms in the 
Ward-Takahashi identity were measured, allowing a determination of the A J5 term. Figure 
I5| shows (ipip), Xtt and AJ5 for a variety of values of L^. Fitting AJ5 to an exponential form 



for Lg = 16 to 40 gives the solid line in the figure and the result 

AJ5 = 0.0096(2) exp(-0.0191(9)L,) xV^of = 6.4/2 (38) 

We see that our data is consistent with AJ5 vanishing as 00, although the decay 

constant is quite small, ~ 1/50. 

Pion poles should dominate the Ward-Takahashi identity when the pions are light and 
the pions should become massless when ruf + m^-es = 0. (This is only strictly true in the 
infinite volume limit.) Thus we look for the pseudoscalar susceptibility in large volumes for 
small total quark mass to behave as 

Xn = a_i/(m/ + mres) + ao + 0{mf + m^es)- (39) 

where the are independent of Ls and mj. This gives a pion pole (for large volumes) at 
rrif = — mres, while oq gives the contribution to the susceptibility of modes whose mass is 
non-zero when the quark mass vanishes. Like {-iptp), ao receives contributions diverging as 
1/a^ and hence may be sensitive to unphysical 5-dimensional modes. For this expression to 
be useful, we do not require the pole term to dominate the remaining terms, but it must 
make a large enough contribution to be visible. 

The A J5 term in Eq. ^ also has a pole contribution coming from the propagation of the 
conventional light pseudoscalar along the s = and Lg — 1 boundaries from to x. This 
light state has non-zero overlap with the midpoint pseudoscalar density for finite Lg, but 
this overlap should be exponentially suppressed. Therefore we expect AJ5 to also have a 
pole at nif = —m^es, giving AJ5 the same form as Xn, namely 

(40) 

Considering the case where the pole terms dominate gives 

(^^) = . L (41 

m + mres 

For (ijjtp) to be finite in this case requires 

a_im + b'_-^ = a_i(m + mrcs) (42) 
so the most general form for A J5 is 

AJ5 = mj-esXTT + bo + 0{mf + rrires), (43) 
Where bo = b'^ — mresQo- Using this then gives 

{iplp) = {nif + mres)X7r + &0 (44) 

up to terms linear in the quark mass. 
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Our procedure for extracting mres from these small volumes involves measuring values for 
X-K and {ipip) for a variety of valence quark masses for a simulation with a fixed dynamical 
quark mass. Since the Ward-Takahashi identity is a consequence of the form of the domain 
wall fermion operator, independent of the weight used to generate the gauge field ensemble 
in which the fermionic observables are measured, it is satisfied by observables measured 
with valence masses. Of course, extrapolations in valence quark mass can lead to problems 
due to the gauge field ensemble including configurations with small fermion eigenvalues that 
are not present when a dynamical extrapolation is done. Here we have a small dynamical 
quark mass present in the generation of the gauge fields, so such effects are expected to be 
unimportant. 



For a given Lg, we simultaneously fit Xvr and {il^i)) to the forms in Eqs. ^ and^. These 
are four parameter fits for oq, a_i,6o and mres and the resulting value for mres we refer to 
as m^^'^^^^. (All measurements of the residual mass from low energy physics should agree. 
We use this notation to detail the explicit technique we have used for this determination.) 
We have used quark masses of 0.02, 0.06, 0.10 and 0.14 in our fits. These fits do not include 
possible correlations between the quantities computed for different values of m/ because the 
correlation matrix itself is poorly determined. 

The results are given in Table pO] , where the errors are all from application of the jack 
knife method. Notice that 6o is negative for all values of L^, meaning that the non-pole 
contributions to A J5 are smaller than mresflo- We have then fit these values of m^^'^'^^^ and 
—60 to the form cq + Ci exp(— aL.;) and found 

- 60 = 0.0104(4) exp(-0.016(2)L,) xV^of = 0.34(19) (45) 
^GMOR) ^ 0.185(6) exp(-0.0280(15)L,) xV^of = 0.28(25) (46) 

(47) 



Figure ITBI shows these values and the fits. 



We can see that both m^^'^*"'^^ and 60 are falling exponentially, but with a very small 
decay constant ~ 1/50. This is in sharp contrast to the decay constant for which is 

1/10. This is further evidence for the distinction between the residual mass that enters 
in low-energy observables and the residual mixing which effects observables dependent on 
degrees of freedom at the cutoff scale. 

Since our determination of the residual mass has been done for small volumes, one can 
worry about the finite volume effects. We have done a similar extraction of the residual 
mass and compared it with determinations of the residual mass from extrapolations of m^ 
for much larger volumes and find reasonable agreement WM. We are continuing to study 
various determinations of the residual mass. 



VI. CONCLUSIONS 

In this work the properties of domain wall fermions relevant to numerical simulations 
of full Nf = 2 QCD at finite temperature were investigated on relatively small lattices of 
size 8^ X 4. Conventional numerical algorithms (the Hybrid Monte Carlo and the conjugate 
gradient) worked without any difficulty beyond the additional computational load of the 
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fifth dimension. Evidence for both confined and deconfined phases was found and the Lg 
and mo dependence of each phase was investigated. 

The domain wall fermion action is expected to preserve the full chiral symmetries of QCD 
for large Lg. For the stronger couplings used for the confined phase simulations, the chiral 
condensate approached its asymptotic value for ^ 32 — 40. However, our determination 
of the residual mass effects present in low energy observables show a residual mass of ~ 0.06 
for Lg = 40. For the weaker couplings needed to study the deconfined, chirally restored 
phase, the residual mass effects are expected to be much smaller for the same Lg, although 
we have not yet measured the residual mass in this region. 

In particular, it was found that for the two flavor theory there is a phase where the 
SU{2) X SU{2) chiral symmetry is broken spontaneously to a full SU{2) flavor symmetry 
and a phase where the full SU{2) x SU{2) chiral symmetry is intact. For the values of Lg we 
used, the dependence of observables on the coupling in the transition region is likely quite 
influenced by the change in the residual mass with the coupling. To suppress this effect 
will require larger values for Lg, thermodynamics studies at larger A^^ (and hence weaker 
coupling) or improved variants of domain wall fermions. 

Our simulations show that domain wall fermions have passed one vital test for numerical 
work, light chiral modes exist at quite strong coupling. A second important result, which was 
expected from work with dynamical fermions in the Schwinger model [^, is that domain 
wall fermions do not present any problems to conventional dynamical fermion numerical 
algorithms. Given these results, we are pursuing simulations of the phase transition on 
larger lattices to achieve more physically meaningful results. The slow falloff of the residual 
mass with Lg can be overcome with more computing power or, hopefully, improvements to 
the formulation. At present, this is all that stands in the way of simulating the Nf = 2 QCD 
phase transition with three degenerate light pions at flnite lattice spacing. 
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APPENDIX A: GAMMA MATRICES 

The Dirac gamma matrices used in this work are: 





1 















( 











-1\ 












i 















1 





7i = 







—i 








, 72 = 







1 










[ 


—i 












\ 


-1 








0/ 




1 








i 


^\ 




( 








1 


°\ 















—i 















1 


73 = 




—i 











, 74 = 




1 


















i 









\ 





1 








18 



75 



/ 1 




V 





1 









-1 








-1/ 



(Al) 



APPENDIX B: EVOLUTION ALGORITHM 

As described in Section ||, we use the Hybrid Monte Carlo algorithm of Gottlieb et al. 
46| extended to include the Pauli-Villars regulator fields. Further, we use a preconditioned 



variant of the Dirac operator specified in Eq. |^ In this Appendix we describe the 



resulting algorithm we use to evolve the gauge fields including the effects of the two flavors 
of domain wall quarks and the Pauli-Villars regulator flelds. 

Following this approach, we generate a Markov chain of gauge flelds Ufj,{x), pseudo- 
fermion flelds ^p, Pauli-Villars flelds $pv and conjugate momenta H^[x) according to the 
distribution: 

Z = j[dU][dH][d^^p][d^F\[d^^Pv][d'^Pv]e''^ (Bl) 

where 

'>^ = Sg + \Y.Ht.{xf +^^p[D^pDf\-^^f +^^py[D^pDFU^.=^^py. (B2) 

Here, the flelds $i? and ^py as well as the preconditioned operator Dp are deflned only on 
odd sites with 

Dp = {b- mo)' - {Dp)oe{Dp)eo (B3) 

where {Dp)oe and {Dp)eo represent the DWF operator of Eq. |] evaluated between odd 
and even or even and odd sites respectively. Note, even and odd are deflned in a flve- 
dimensional sense, e.g. for an even site the sum of all flve coordinates is an even number. 
Eq. employs the usual preconditioning scheme for Wilson fermions implemented in 5 



dimensions. Similar considerations justify the form used for the Pauli-Villars action. Since 
detDp = det{(5 — mo)Dp}, we have rescaled both the flelds $i? and to introduce the 
extra factor of (5 — mo) into Eq. ^3| in order to simplify the subsequent algebra. 

To begin a new HMC trajectory, we start with the values of the gauge flelds U^{x) pro- 
duced by the previous trajectory. We then choose Gaussian distributed flelds s)i?, 
rj{x,s)pv and H^{x) from which we construct the flelds ^p = Dprjp and ^pv = 
{Dp^\mf=i)r]pv- Here we have introduced new fleld variables H^{x), conjugate to the link 
matrices, which are elements of the algebra of SU{3), and hence traceless and hermitian. 

Next, we carry out the molecular dynamics time evolution of the flelds H^{x) and U^{x) 
according to equations of motion which are phase space volume preserving and conserve the 
flctitious 6-dimensional "energy" Ti of Eq. [B2| . The flrst of these Hamilton-like equations 
determines the relation between U^{x) and the conjugate variable H^[x): 
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dr 



(B4) 



The second equation can be derived from the requirement that Ti. is r-independent. First, 
following Gottlieb et al. PB| one writes: 



d7 



ETr 



iH,{x)F,{x) + -^^H,{x) 

(XT 



(B5) 



Then the constancy of Ti is insured if for the second equation of motion we impose: 

The subscript TA indicates the traceless anti-hermitian part of the matrix, a restriction 
required by the traceless, hermitian character of the variables H^{x). (The definition of 
F^{x) implied by Eq. ^5| makes F anti-hermitian and it is only the traceless part of F that 
enters that equation.) 

Finally we will determine the specific form for the force term F^{x). This can be done 
by using the general formula 



d_ 
d7 



{^'\Df\ 



'(x, s)^H^{x)U,{x){l T in^{x + /i, s) 

x,s 

-ij'{x + /i, s)^U^{x)^H^{x){l ± 7^)V^(a;, s)] 



(B7) 



which follows immediately from Eq's. ^ and where the lower choice of signs corresponds 
to the case olD^. Now we re-express the derivative: 



d_ 
d^' 



Xf, 



(B8) 



where we construct $ = Dprjp from the Gaussian source rjp and then obtain xf by solving 
DpDpXF = ^F- Now we must evaluate 



x'f 



d 

d^ 



^D^Dp 



Xf = ^{xf\ [(5 - mo)' - (/^^)oe(/^^)eo] 
■[{5-moy-{Dp)UDp),o]\XF), 



(B9) 



We will obtain eight terms by letting the derivative act on each of the four Dp operators. 
Four of those terms will involve Ufj_{x) and four U^^xY , with the final four terms being the 
hermitian conjugates of the first four. Combining Eq.'s |B7| , P8| and [B9| , we find: 

tT{H^{x)F^{x)} = W {xf{x, s)H^{x)U,{x){1 + ^^{{x + /x, s)\{D^p)eoDp\xF) 

^ s 

+ {xF\{D^p)oe\x, s)H^{x)U^{x){l + 7'^)(X + ^Ji,s\Dp\xF) 
+ {xf\D^p\x,s)H^{x)U^{x)(1--i^){x + li,s\{Dp)eo\XF) 
+ {xF\D^F{Dp)oe\x,s)H^{x)U^{x){l-^^')x{x + ^l,s)-\i.c.]. (BIO) 
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This expression can be written in a very simple form if we define two new spinor quantities: 

^ I -{x,s\{D^F)eoDF\XF) {x,s) even ^^^^^ 
\ -{x,s\Df\xf) {x,s) odd 

^('3, {x,s\{DF)eo\XF) {x,s) even ,^^2) 

I XFix,s) ix,s) odd 

Using these quantities in Eq. |B10| and factoring out the generator H^{x) gives: 

F[F\^,{x) = -^f//.(a;)^trspin Ul - j^,)v{x + p., s)w\x, s) 

+ (1 + 7^)w(x + /i, s)f ^(x, s) — h.c. (B13) 

where we have added now the subscript [F] to distinguish this fermion force from that 
produced by the Pauh-Villars fields described below. Since there are no gauge fields in the 
extra direction, it is not surprising that this looks very similar to the Wilson fermion force 
with an additional sum over the s-direction. 

The force term produced by the Pauli-Villars fields is closely related to that derived 
above. We need only replace the field xf with ^py, set mj = 1 and change the sign of the 
resulting force: 

F,FvU-) = -FiFU-)i^^,,,^.^^^- (B14) 
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FIG. 2. Valence extrapolations of (V'V') for = 0-1) = 1-65 and Ls = 8. The circles are 
for /? = 5.0, the squares for /? = 5.4. 
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FIG. 3. The simulation "time history" of {tpi') for /? = 5.0 (upper graph) and /3 = 5.4 lower 
^raph for m/ = 0.1, mo = 1.90 and Lg = 8. 
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FIG. 4. Valence extrapolations of (V'V') ^or mj = 0.1, mo = 1.90 and Ls = 8. The circles are 
results for P = 5.0, the squares are for /? = 5.4. 
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FIG. 5. Full QCD extrapolations of (^V') for mo = 1.90 and Lg = 8. The circles are results 
for (3 = 5.0, the squares are for /3 = 5.4. The dashed line is the fit to the quenched data given in 
Figure |^ while the solid line is a fit to the dynamical values for . 
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FIG. 7. The simulation "time history" of (i/^V) ^^r (3 = 5.2 (upper graph) and [3 = 5.45 (lower 
graph) for mf = 0.02, mo = 1.9 and Lg = 16. The initial configuration was chosen in the opposite 
phase, i.e. ordered for (3 = 5.2 and disordered for (5 = 5.45. 
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FIG. 9. Full QCD values for {t/ji/j) for mo = 1.90 and P = 5.2 plotted versus Lg for different 
values oiruf. The curves are fits of the form cq + ci exp(— aL^). 
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FIG. 10. Full QCD values for (V'V') for mo = 1.90 and /3 = 5.45 versus nif. The circles are for 
Lg = 8, the squares for 16. 
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FIG. 11. Full QCD values for (V'V') for mo = 1.90 and /? = 5.45 plotted versus Lg for different 
values of ruf. The curves are fits of the form cq + ci exp(— aL^). 
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FIG. 13. Full QCD results for (i/'^) for = 0.1 and Lg = 12 for different values of mo and /3 

near the transition region. Note that the vertical scale decreases by a factor of five as tuq decreases 
from 2.4 to 1.15. This is needed to follow the large decrease in the scale of (V'V') which results from 
a combination of the decreasing lattice spacing that follows from increasing /3 and the diminishing 
overlap of the light fermion states with the walls. 
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FIG. 14. The critical value of /3 as a function of mo- The line is the value for a 24^ x 4 lattice. 
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FIG. 15. The chiral condensate, {il^ip), the pion susceptibiHty, Xtt s^nd AJ5 versus Lg for (3 = 5.2, 
mo = 1.9 and = 0.02. 
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TABLES 





Ls 


tra^ 


. len. 


# traj. 


acc. 


/ AFf \ 


(plaq) 


{\w\) 


/-r , \ 

(V'V') 


0.02 


8 


1 

64 


X 32 


400-800 


0.89 


0.98(1) 


0.456(2) 


0.061(2) 


0.0133(2) 




10 


1 

64 


X 32 


200-2000 


0.86 


0.995(8) 


0.4460(9) 


0.048(2) 


0.0124(1) 




12 


1 

64 


X 32 


200-2000 


0.84 


1.01(1) 


0.4428(6) 


0.048(2) 


0.01123(7) 




16 


1 

64 


X 32 


550-2000 


0.75 


0.98(2) 


0.4388(9) 


0.049(2) 


0.00987(9) 




24 


1 
100 


X 50 


350-2000 


0.73 


0.95(2) 


0.4359(7) 


0.047(3) 


0.0088(1) 




32 


1 
100 


X 50 


300-2000 


0.72 


1.03(2) 


0.4317(7) 


0.045(2) 


0.00835(7) 




40 


1 

128 


X 64 


300-1350 


0.73 


1.02(3) 


0.4342(6) 


0.044(2) 


0.00772(8) 


0.06 


8 


1 

50 


X 25 


200-950 


0.83 


0.98(1) 


0.450(1) 


0.046(3) 


0.0176(2) 




1(3 


1 

f)l 


X 32 


200-820 


0.84 


0.99(2) 


0.1361(8) 


0.015(3) 


0.0135(1) 


0.1 


8 


1 

40 


X 20 


300-800 


0.57 


0.89(5) 


0.4437(8) 


0.040(3) 


0.02109(7) 




10 


1 

50 


X 25 


200-800 


0.83 


1.00(3) 


0.4405(6) 


0.036(2) 


0.01927(9) 




12 


1 

40 


X 20 


400-800 


0.43 


1.2(3) 


0.437(1) 


0.032(2) 


0.01838(4) 




16 


1 

64 


X 32 


200-800 


0.80 


0.98(2) 


0.435(2) 


0.035(2) 


0.01709(9) 




24 


1 

64 


X 32 


200-800 


0.72 


0.94(2) 


0.433(1) 


0.033(2) 


0.01596(7) 




32 


1 
100 


X 50 


200-800 


0.82 


0.99(2) 


0.4305(5) 


0.037(2) 


0.01547(7) 




40 


1 
100 


X 50 


200-800 


0.78 


1.00(5) 


0.432(1) 


0.035(2) 


0.01524(5) 


0.14 


8 


1 

40 


X 20 


200-860 


0.63 


1.03(6) 


0.4433(7) 


0.033(5) 


0.0241(1) 




16 


1 

64 


X 32 


200-800 


0.85 


1.03(2) 


0.433(1) 


0.030(1) 


0.02017(6) 


0.18 


8 


1 

40 


X 20 


200-1200 


0.70 


1.02(3) 


0.4410(7) 


0.030(1) 


0.02686(7) 




If) 


1 

()1 


X 32 


200-800 


0.81 


0.98(2) 


0.-132(1) 


0.033(1) 


0.02309(5) 



TABLE I. Data for an 8^ x 4 lattice with (3 = 5.2, mo = 1.9 



40 



ruf 


Ls 


tra 


. len. 


# traj. 


acc. 


/ -AH"\ 


(plaq) 


{\w\) 


/-r 1 \ 


0.02 


8 


1 

64 


X 32 


200-800 


0.91 


1.005(7) 


0.5376(7) 


0.226(4) 


0.00415(6) 




10 


1 

64 


X 32 


200-1000 


0.91 


0.992(8) 


0.5328(6) 


0.207(4) 


0.00319(5) 




12 


1 

64 


X 32 


200-800 


0.95 


1.009(9) 


0.5300(4) 


0.202(5) 


0.00270(3) 




16 


1 

64 


X 32 


200-800 


0.90 


1.02(1) 


0.5266(8) 


0.199(4) 


0.00237(6) 




24 


1 

64 


X 32 


400-1200 


0.86 


0.98(2) 


0.5257(7) 


0.187(3) 


0.00216(6) 




32 


1 
100 


X 50 


400-800 


0.94 


1.00(2) 


0.524(2) 


0.180(5) 


0.00209(5) 


0.06 


8 


1 

50 


X 25 


200-1000 


0.86 


0.99(3) 


0.536(1) 


0.217(3) 


0.0080(1) 




10 


1 

64 


X 32 


200-1000 


0.92 


0.994(7) 


0.5313(6) 


0.203(4) 


0.00704(5) 




12 


1 

64 


X 32 


200-1000 


0.89 


1.013(8) 


0.5286(8) 


0.195(4) 


0.00666(5) 




16 


1 

64 


X 32 


400-800 


0.76 


1.02(4) 


0.525(2) 


0.192(4) 


0.00637(7) 




24 


1 

64 


X 32 


300-1000 


0.84 


1.00(1) 


0.521(2) 


0.174(6) 


0.00617(9) 




32 


1 

64 


X 32 


500-1000 


0.80 


1.00(2) 


0.525(2) 


0.189(3) 


0.00592(4) 


0.1 


8 


1 

50 


X 25 


300-800 


0.83 


0.98(2) 


0.5336(6) 


0.211(4) 


0.01174(4) 




10 


1 

50 


X 25 


300-990 


0.88 


0.99(1) 


0.5310(9) 


0.200(2) 


0.01075(5) 




12 


1 

50 


X 25 


600-1200 


0.74 


1.01(4) 


0.528(1) 


0.197(4) 


0.01838(4) 




16 


1 

64 


X 32 


400-800 


0.79 


1.01(3) 


0.523(1) 


0.170(5) 


0.0103(1) 




24 


64 


X 32 


400-2000 


0.86 


0.991(8) 


0.512(1) 


0.170(8) 


0.0102(1) 




32 


1 

64 


X 32 


300-1000 


0.81 


0.98(2) 


0.519(1) 


0.159(5) 


0.01011(9) 


0.14 


8 


1 

50 


X 25 


200-800 


0.83 


1.01(1) 


0.533(1) 


0.210(3) 


0.01531(9) 




16 


1 

64 


X 32 


600-1200 


0.76 


0.98(2) 


0.520(1) 


0.159(9) 


0.0143(1) 


0.18 


8 


1 

50 


X 25 


400-800 


0.81 


1.03(2) 


0.5314(6) 


0.202(4) 


0.01884(5) 




16 


1 

64 


X 32 


600-1200 


0.78 


0.94(2) 


0.515(1) 


0.141(8) 


0.0182(2) 



TABLE II. Data for an 8^ x 4 lattice with (3 = 5.45, mo = 1.9 
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HMC traj 


. len: 




CG stop cond: 


10-6 




P 


start 


# traj 


acc. 


/„-ah\ 
\^ / 


(plaq) 




\VV 1 


5.45 





100-800 


0.87 


0.99(1) 


0.470(1) 


0.0168(6) 


0.00276(1) 


5.55 





200-800 


0.87 


0.98(1) 


0.4933(6) 


0.023(1) 


0.002916(6) 


5.65 





300-800 


0.87 


1.00(2) 


0.5218(9) 


0.054(6) 


0.00305(2) 


5.75 


D 


300-800 


0.86 


0.986(9) 


0.5571(7) 


0.196(7) 


0.002875(7) 


5.85 


D 


300-800 


0.85 


1.00(1) 


0.5719(7) 


0.234(3) 


0.002881(3) 


5.95 


D 


200-800 


0.87 


0.99(1) 


0.5857(5) 


0.262(2) 


0.002898(3) 



TABLE III. Data for an 8^ x 4 lattice with mo = 1.15, Lg = 12, and m/ = 0.1. 



HMC traj. len: ^ x 25, CG stop cond: 10-' 



/? 


start 


# traj 


acc. 




(plaq) 


(iw^l) 




5.35 





100-800 


0.87 


1.01(1) 


0.4435(7) 


0.0189(8) 


0.00497(1) 


5.45 





100-800 


0.86 


0.99(1) 


0.4630(7) 


0.0215(8) 


0.00522(1) 


5.55 





300-800 


0.86 


1.00(1) 


0.487(1) 


0.032(3) 


0.00539(3) 


5.65 


D 


400-800 


0.86 


1.00(2) 


0.540(2) 


0.180(6) 


0.00457(4) 


5.75 


D 


300-800 


0.85 


0.99(1) 


0.5598(8) 


0.224(3) 


0.00445(1) 


5.85 


D 


200-800 


0.89 


1.011(8) 


0.5744(3) 


0.254(3) 


0.004409(5) 




TABLE IV. Data for an 8^ x 4 lattice with 


mo = 1.4, Ls = 


12, and ruf = 


0.1. 






HMC traj 


• len: ^ 


X 25, 


CG stop cond: 10-^ 






start 


# traj 


acc. 


(e— ) 


(plaq) 


{\w\) 




5.25 





200-800 


0.82 


0.97(2) 


0.4289(5) 


0.027(2) 


0.01000(2) 


5.35 





400-800 


0.68 


0.98(4) 


0.451(3) 


0.035(4) 


0.01000(9) 


5.45 


D 


400-800 


0.74 


1.09(5) 


0.4769(8) 


0.049(7) 


0.00985(7) 


5.55 


D 


600-1200 


0.80 


1.00(3) 


0.531(1) 


0.175(7) 


0.00718(7) 


5.65 


D 


400-800 


0.79 


0.96(2) 


0.5507(9) 


0.216(4) 


0.00677(2) 


5.75 


D 


200-800 


0.88 


0.989(7) 


0.5663(4) 


0.243(3) 


0.00658(1) 



TABLE V. Data for an 8^ x 4 lattice with mo = 1.65, Lg = 12, and m/ = 0.1. 
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HMC traj. len: ^ x 25, CG stop cond: 10"^ 

/? start # traj acc. {^~^") (plaq) {\W\) (j)^') 

5I5 O 200-800 083 1.02(2) 0.4191(8) 0.029(1) 0.01485(5) 

5.25 O 400-800 0.66 0.97(5) 0.4381(6) 0.038(2) 0.01458(5) 

5.35 O 400-800 0.63 0.97(5) 0.471(2) 0.052(3) 0.0134(2) 

5.45 O 400-800 0.76 1.01(3) 0.515(2) 0.161(4) 0.0097(1) 

5.55 D 400-800 0.79 1.05(5) 0.540(1) 0.200(9) 0.0088(1) 

5.65 D 200-800 0.89 1.01(2) 0.5570(5) 0.242(4) 0.00828(2) 

TABLE VI. Data for an 8^ x 4 lattice with mo = 1.8, Lg = 12, and m/ = 0.1. 

CG stop cond: IQ-^ 

13 start traj. len. # traj acc. {^~^") (P^^l) l^'^) 

YX) O ^ X 20 200-800 037 0.8(1) 0.4002(8) 0.032(2) 0.01919(5) 

5.2 O ^ X 20 400-800 0.43 1.2(3) 0.437(1) 0.032(2) 0.01838(4) 

5.25 O ^ X 25 400-800 0.65 1.10(9) 0.452(1) 0.049(6) 0.0174(2) 

5.35 D 35 X 25 600-1200 0.69 0.95(5) 0.493(2) 0.107(9) 0.0135(4) 

5.45 D 35 25 600-1200 0.74 1.01(4) 0.528(1) 0.197(4) 0.01039(7) 

5.55 D ^ 25 400-830 0.82 1.00(1) 0.5463(5) 0.227(6) 0.00974(4) 

5.65 D in ^ 25 400-800 088 1.03(1) 0.5613(8) 0.248(5) 0.00943(4) 

TABLE VII. Data for an 8^ x 4 lattice with tuq = 1.9, = 12, and m/ = 0.1. 

HMC traj. len: ^ x 25, CG stop cond: 10"^ 

P start # traj acc. {^~^^) (plaq) {\W\) (^V>) 

K05 O 200-800 077 0.99(3) 0.4192(8) 0.035(1) 0.02324(8) 

5.15 O 200-800 0.75 0.98(3) 0.442(1) 0.042(3) 0.0215(2) 

5.25 O 200-1200 0.79 1.03(1) 0.474(1) 0.080(7) 0.0181(3) 

5.35 D 200-800 0.83 1.00(2) 0.5130(7) 0.173(6) 0.0130(2) 

5.45 D 200-800 0.87 1.02(2) 0.5349(5) 0.203(3) 0.01157(4) 

5.55 D 200-800 0.85 1.01(1) 0.5503(4) 0.235(3) 0.01099(2) 

TABLE VIII. Data for an 8^ x 4 lattice with mo = 2.0, Lg = 12, and m/ = 01. 
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HMC traj. 


len: ^ 


X 25, 


CG stop cond: 


10-6 




P 


start 


• 

# traj 


acc. 




(plaq) 


/I Wl\ 




4.85 





200-800 


0.69 


0.95(3) 


0.4004(8) 


0.034(2) 


0.0323(2) 


4.95 





200-800 


0.72 


0.97(3) 


0.419(2) 


0.040(2) 


0.0302(3) 


5. (Jo 


(J 


zUO-oOU 


0.48 


0.92(4] 


0.443(2) 


0.052(3) 


0.0272(5) 


5.15 





200-1200 


0.62 


1.01(3) 


480('3') 


0.12(1) 


0208('7'l 


5.25 





400-800 


0.70 


0.97(5) 


0.5105(4) 


0.185(2) 


0.01559(8) 


5.35 


D 


400-800 


0.69 


1.01(4) 


0.529(1) 


0.216(6) 


0.0141(1) 


5.45 


D 


400-800 


0.71 


1.00(3) 


0.5453(7) 


0.230(4) 


0.01330(5) 




TABLE IX. Data for an 


8^ X 4 lattice with 


mo = 2.15, Ls = 


: 12, and m/ 


= 0.1. 



HMC traj. len: ^ x 25, CG stop cond: lO'^ 





start 


# traj 


acc. 


(e— ) 


(plaq) 


m\) 




4.65 





100-800 


0.63 


1.02(5) 


0.3953(6) 


0.046(3) 


0.0484(3) 


4.75 





200-800 


0.68 


0.99(5) 


0.4156(9) 


0.054(3) 


0.0442(3) 


4.85 





300-800 


0.70 


0.94(4) 


0.439(2) 


0.069(5) 


0.0380(6) 


4.95 





200-800 


0.77 


1.02(4) 


0.4779(6) 


0.155(4) 


0.0257(2) 


5.05 


D 


200-800 


0.80 


1.01(2) 


0.4987(9) 


0.190(2) 


0.0220(2) 


5.15 


D 


200-800 


0.84 


1.01(3) 


0.5170(5) 


0.221(3) 


0.01962(7) 




TABLE X. Data for 


an 8^ X 4 lattice with mo 


= 2.4, Ls = 


12, and ruf = 


0.1. 



Ls mrcs X^/dof 

10 0.149(5) 0.0094(5) 0.8(4) 

12 0.129(2) 0.0080(2) 1.6(4) 

16 0.113(3) 0.0080(4) 1.1(5) 

24 0.095(2) 0.0075(3) 1.5(7) 

32 0.078(2) 0.0068(5) 0.7(4) 

40 0.059(3) 0.0048(3) 1.7(9) 



TABLE XI. Values for mrcs and —bo versus Lg from fits to valence quark data with the 
dynamical quark mass fixed at nif = 0.02. 
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